% Load county-to-county pairs data;
% SCI_data.csv is the county-county-SCI datasets;
clear all
% Directed graphs;
R1 = csvread('SCI_data.csv', 1);
n = length(unique(R1(:, 1)));

% Code county fipes from 1 to n;

% Extract the 'own_county' and 'friend_county' columns
own_county = R1(:, 1);
friend_county = R1(:, 2);

% Create a unique list of all FIPS codes
all_counties = unique([own_county; friend_county]);

% Initialize matrices to store the IDs
own_county_id = zeros(size(own_county));
friend_county_id = zeros(size(friend_county));

% Replace FIPS codes with IDs
for i = 1:length(all_counties)
    own_county_id(own_county == all_counties(i)) = i;
    friend_county_id(friend_county == all_counties(i)) = i;
end

% add the 'own_county' and 'friend_county' columns in the R1 matrix with the new ID columns
R1(:, 4) = own_county_id;
R1(:, 5) = friend_county_id;

% Construct adjacency matrix
ADJ = accumarray([R1(:,4),R1(:,5)],R1(:,3),[n n]);

% ADJ is already symmetric: no, some small scis are rounded off the other way;
sum(sum(ADJ ~= ADJ'));
ADJ_sm = max(ADJ, ADJ');
clear ADJ;
% eigenvector centrality
[D_2016, V_2016] = eig(ADJ_sm); 
dominant = find(diag(V_2016)==max(diag(V_2016)));
eigenvector_2016 = D_2016(:, dominant)/sum(D_2016(:, dominant));
% degree centrality
outdegree_2016 = sum(ADJ_sm, 2) - diag(ADJ_sm,0);
outdegree_2016 = outdegree_2016./sum(outdegree_2016);
% Information centrality; 
B = diag(sum(ADJ_sm,2)) - ADJ_sm + ones(size(ADJ_sm));
B_inverse = B^(-1);
clear B
I = zeros(size(ADJ_sm));
n = length(ADJ_sm(1,:));
for i = 1:n
    for j = 1:n
        I(i,j) = B_inverse(i,i) + B_inverse(j,j)- 2*B_inverse(i,j);
    end;
end;
% IC as harmonic means of I;
% need to change diag(I) to be inf;
IC_2016 = harmmean(I.^(-1),2);
IC_2016 = IC_2016./sum(IC_2016); 

centrality_2016 = [all_counties, IC_2016, outdegree_2016, eigenvector_2016];

%% Add county fips and save;
hdr = 'county, infc, outdegree, eigenvector';
dlmwrite('centrality.csv', hdr, '');
dlmwrite('centrality.csv', centrality_2016, '-append', 'delimiter', ',', 'precision', 20); 
